

rbeta(100, 2, 5)
# Plotting the beta distribution  
x <- seq(0, 1, length = 100)
y <- dbeta(x, 0.2, 5)

plot(x, y, type = "l", lwd = 2, col = "blue",
     main = "Beta Distribution (shape1=2, shape2=5)",
     xlab = "x", ylab = "Density")



beta_area(0.4, 0.6, c(7,10))


beta.select(list(x=0.55, p=0.5),
            list(x=0.8, p=0.9))

x1<-0.55
p1<-0.5
x2<-0.8
p2<-0.9

objective <- function(params){
  alpha <- params[1]
  beta <- params[2]
  
  err1 <- pbeta(x, alpha, beta) - p1
  err2 <- pbeta(x, alpha, beta) - p2
  
  return (err1^2 + err2 ^2)
}

start <- c(2,2)

result <- optim(start, objective, method ="L-BFGS-B",
                lower = c(1e-3, 1e-3))

alpha_est <- result$par[1]
beta_est <- result$par[2]

logit1 <- log(p1 / (1-p1))
logit2 <- log(p2 / (1-p2))

b1 <- (logit2 - logit1) / (x2 - x1)
a1 <- logit1 - b1 * x1

print(b1)
print(a1)




